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1.  Introduction 


Given  a  larjge  set  of  measurements  of  some  quantity  or  variable,  it 
is  often  useful  to  model  the  data  using  some  statistical  distribution 
function.  For  example  if  one  has  records  of  the  "visibility"  at  Milden- 
hall,  England  for  10  a.m.  February  over  a  number  of  years,  one  may  fit 
a  Weibull  distribution  to  the  data.  The  Weibull  distribution  has  two 
Parameters,  and  the  values  selected  for  the  two  parameters  are  the  ones 
for  which  the  model  best  fits  the  data. 

This  report  documents  a  FORTRAN  program  that  has  been  written  to 
estimate  the  parameters  of  a  statistical  distribution  function  which 
best  fits  a  set  of  measurements  on  some  variable.  The  fit  is  "best" 
in  the  sense  that  the  model  cumulative  distribution  function  and  the 
empirical  cumulative  distribution  function  (from  the  data)  are  closest 
to  each  other  in  the  least  squares  sense. 

Suppose  the  measurements  are  ordered  from  smallest  to  largest, 
that  is  x^  <  x^2j  ^  x(3)  x^  where  x^  represents  the  ith 

smallest  measurement.  Then  the  empirical  cumulative  distribution  function 
may  be  defined  as 

Vx)  “  for  x(i>  - x "  xu+i) 

«  0  for  x  <  x^  .  (1.1) 


If  F(x;9)  is  the  model  cumulative  distribution  function,  then  the  values 
for  9  (9  may  be  a  vector  of  values)  which  are  chosen  are  those  which 
minimize  the  expression 

N  r  m2 

l  [(2i-l)/(2S)-F(xi;Q)j  .  (1.2) 

In  the  FORTRAN  program,  the  determination  of  9,  is  accomplished  by  non¬ 
linear  regression  techniques  where  FN(x)  =  (2i-l)/(2N)  for  i  »  1,2,...,N, 

is  the  dependent  variable.  Expression  (1.2),  the  quantity  to  be  minimized, 
is  the  "Residual  Sum  of  Squares".  A  detailed  description  of  the 
techniques  used  to  fit  distributions  to  data  using  non-linear  regression 
techniques  is  given  in  Heuser,  Somerville  and  Bean  (19803 
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2.  Flow  Chart 


In  non-linear  regression,  the  model  may  be  written 

y  -  F(x;0)  +  e  ,  (2.1) 

where  3  represents  a  vector  of  unknown  parameters. 

The  usual  technique  is  to  linearize  F(x;0)  by  the  use  of  a  first  order 
Taylor  Series  expansion  about  guessed  values  0O.  The  expression  (2.1) 
is  then  linear  in  (0-0Q),  and  the  usual  least  squares  regression 
methods  may  be  used  to  estimate  0-0Q,  the  "correction"  to  the  original 
guessed  value.  The  procedure  is  then  repeated  with  a  first  order 
Taylor  Series  expansion  about  the  "corrected"  guessed  value  for  0, 
the  process  terminating  when  the  percentage  reduction  in  residual  sum 
of  squares  is  less  than  some  specified  amounts.  The  flow  chart  below 
outlines  the  program. 
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3.  FORTRAN  Code 


C***************************** ****************************** *********** 

c  * 

c  title:  weibull  nonlinear  regression  program  # 

C  THE  FOLLOWING  PROGRAM  REGRESSES  VISIBILITY  DATA  ON  THE  * 

C  WEIBULL  DISTRIBUTION.  * 

C  * 

C  INPUT:  N.  X'liN),  ACC  ,  AND  Y  < 1 J  N  >  <SEE  VARIABLE  DICTIONARY)  * 

C  N,  X'lJN),  AND  ACC  ARE  INPUT  ONCE  IN  THE  BEGINNING  OF  THE  * 

C  PROGRAM.  Y(1:n>  IS  INPUT  ONCE  FOR  EACH  REGRESSION.  * 

C  * 

C  OUTPUT:  SUMMARY  STATISTICS  OF  EACH  REGRESSION  INCLUDING  ALPHA,  BETA,* 
C  MONTH,  HOUR,  SID,  RMS,  COUNT,  X(1JN),  Y(1JN),  PREDICTED  * 

C  VALUES  OF  THE  DISTRIBUTION,  AND  THE  RESIDUALS.  * 

C  * 

c  subroutines:  guess,  sse,  weibull,  pssea,  psseb,  SUCCESS,  FAIL,  * 

c  SECANT,  CORVFC ,  MODL.IN  * 

C  * 

C  FOP  A  GENERAL  OVERVIEW  OF  THE  REGRESSION  PROBLEM,  SEE  ’LEAST  * 

SQUARES  FITTING  OF  DISTRIBUTIONS  USING  NON-LINEAR  REGRESSION’  BY  * 
MARK  HEUSER,  PAUL  SOMERVILLE,  AND  STEVE  BEAN.  * 

* 

******** *******************************************  ************* ****** 

********************************************************************** 

* 

VARIABLE  DICTIONARY  * 

* 

NJ  THE  NUMBER  OF  OBSERVATIONS  <MAXIMUM  OF  15)  * 

X  < 1 : N  >  :  THE  VALUES  OF  THE  INDEPENDENT  VARIABLE  (DISTANCE)  * 

c  yu:n>:  the  observed  visibility  probabilities  at  each  x  * 

C  ALPHA, BETA:  THE  PARAMETERS  IN  THE  WEIBULL  MODEL  * 

C  3TARTA, STARTS:  THE  STARTING  VALUES  FOR  ALPHA  AND  BETA  COMPUTED  BY  * 
C  THE  SUBROUTINE  'GUESS'  * 

C  CORA,CORB:  THE  CORRECTION  VECTORS  FOR  ALPHA  AND  BETA  COMPUTED  BY  * 

C  THE  SUBROUTINE  'CORVEC'  * 

C  NRSS,ORSS:  THE  RSS  FOR  TWO  CONSECUTIVE  ESTIMATES  OF  ALPHA  AND  * 

C  BETA.  NRSS  IS  FROM  THE  NEWER  ESTIMATE?  ORSS  I.S  FROM  * 

C  THE  OLDER  ESTIMATE.  * 

C  MONTH, HOUR, SID?  MONTH,  HOUR,  AND  STATION  IDENTIFIERS  * 

c:  count:  a  loop  counter  * 

C  CONVERGE:  a  LOGICAL  VARIABLE  THAT  INDICATES  CONVERGENCE. 

C  ACC:  AM  INTEGER  VALUE  CONTROLLING  THE  ACCURACY  OF  THE  STARTING 
':  VALUES  FOR  ALPHA  AND  BETA.  SEE  SUBROUTINE  'GUESS'. 

C 

r.  *************************************************  ******************** 
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c 

c 


c 

10 


REAL  ALPHA  * BETA  *  X ( 15 ) » Y( 15) *  START A * STARTS . NRSS . ORSS * CORA * CORB 

REAL  SSE  t  WE  I BUL  *  RMS 

INTEGER  S ID. MONTH* HOUR. N* COUNT* ACC 

LOGICAL  CONVERGE 

COMMON  N.X.Y 

l-J  R I T  E  (  6  *  2  0  0  )  !  PRINT  A  TITLE 


READ  *  N 

READ  f <X( I ) * 1  =  1 »N> 
READ » ACC 


INPUT  THE  NUMBER  OF  OBSERVATIONS 

INPUT  VALUES  Of-  THE  INDEPENDENT  VARIABLE 

INPUT  LEVEL  OF  ACCURACY  OF  STARTING  VALUE 


1  THE  FOLLOWING  LOOP  INPUTS  AND  REGRESSES  ON  THE  EMPIRICAL 
1  DISTRIBUTION.  THE  LOOP  (AND  THE  PROGRAM;  TERMINATES  ON 
1  END  OF  FILE. 


READ  (5*100*  END  =40 )  (Y(I>.I=1*N) *  S I D  *  MONTH  *  HOUR 


CALL  GUESS(STARTA.STARTB*ACC>  )  GET  STARTING  VALUES  FOR  ALPHA 
ALPHA=STARTA  !  AND  BETA 

BETA=STARTB 
NR3S=SSE( ALPHA * BETA) 


!  THE  FOLLOWING  LOOP  SOLVES  FOR  ALPHA  AND  BETA.  CONVERGENCE  IS 
1  ASSUMED  WHEN  THE  PROPORTIONAL  CHANGE  IN  THE  RSS  FOR  TWO  CON 
!  SECUTIVE  ESTIMATES  IS  LESS  THAN  IE-7. 

C 

C0UNT=0  !  INITIALIZE  THE  LOOP  CONTROL 

CONVERGE-. FALSE.  !  VARIABLES 

C 

20  IF  < (COUNT. GT. 50) .OR. (CONVERGE) )  GOTO  30 

ORSS-NRSS 

CALL  CORVEC ( ALPHA  *  BETA  * CORA  *  CORD  > 

CALL  MODLIN( ALPHA * BETA f CORA .CORK) 

NRSS=SSE( ALPHA . BETA) 

CONVERGE~ABS( ORSS-NRSS) ,LT. ( ORSS* 1 . OE-7 > 

C0UNT=C0UNT+1 
GOTO  20 
r 

30  IF  (CONVERGE)  THEN 

CALL  SUCCESS (SID* MONTH  *  HOUR  *  ALPHA * BETA  *  NRSS *  COUNT > 

ELSE 

CALL  FAIL (SID* MONTH  *  HOUR* ST ART A, STARTS .ALPHA* BET A* NRSS  ) 

END  IF 

r 

GOTO  10 
C 

AO  STOP 

r 

i  •.'•0  FORMAT  (IX  .  1  4FA  .3*15*12*11) 

COO  F0PMAT(///*35X* 'NONLINEAR  REGRESSION  OF  THE  UEIBULL  MODEL  ON 
%  'VISIBILITY  DATA'  .///) 

END 
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c  ■ 


C  # 
C  WEIBUL  IS  A  REAL  FUNCTION  THAT  COMPUTES  THE  VALUE  OF  THE  UEIBULL  * 
C  DISTRIBUTION  FOR  THE  INPUT  PARAMETERS  X,  ALPHA  ,  AND  BETA.  ALL  * 
C  COMMUNICATION  WITH  THE  PROCEDURE  IS  THROUGH  THE  PARAMETER  LIST  # 
C  AND  FUNCTION  NAME.  * 
C  * 


C  *  *  *  **  *  ***  ##  *  *  DC  ************  *  *##  *  ###*#***#**  ******  *  *  #***##*#*##**  **  *##  * 

C 

REAL  FUNCTION  WE IBUL ( X , ALPHA » BET A ) 

REAL  X  , ALPHA , BETA 

WEIBUL  =  1 . O-EXP (-ALPHA* <X**BET A)  ) 

RETURN 

END 

C 

c *****************  **************************************************** 

c  * 

C  SSE  IS  A  REAL  FUNCTION  THAT  COMPUTES  THE  SUM  OF  THE  SQUARED  ERRORS  * 
C  IN  THE  UEIBULL  MODEL  AS  A  FUNCTION  OF  ALPHA  AND  BETA.  COMMUNICA-  * 
TION  WITH  THE  PROCEDURE  IS  DONE  THROUGH  THE  PARAMETER  LIST,  THE  * 
FUNCTION  NAME?  AND  THE  COMMON  BLOCK.  * 

* 

C*  **********************  *#***#**#*#**#********#)|[  ********************** 

C 

REAL  FUNCTION  SSE ( ALPHA , BET A ) 

INTEGER  I , N 

REAL  ALPHA, BET A, X< 151 , Y<15> , WE IBUL 
COMMON  N  ?  X  ,  Y 
SSE=0.0 
DO  10  1=1  ,  N 

SSE -SSE+ <  Y ( I > -WEIBUL  <  X ( I ) , ALPHA , BET A ) >**2 
1.0  CONTINUE 

RETURN 
END 

r 

C ******* ******************************* ******* ****************** ****** 


c  CORVEC  IS  A  SUBROUTINE  THAT  COMPUTES  THE  CORRECTION  VECTORS  CORA  * 
C  AND  CORB.  COMMUNICATION  WITH  THE  PROCEDURE  IS  DONE  THROUGH  THE  * 
C  PARAMETER  LIST  AND  THE  COMMON  BLOCK.  THE  INPUT  PARAMETERS  ARE  * 
C  ALPHA  AND  BET  A  j  OUTPUT  PARAMETERS  ARE  CORA  AND  CORB,  * 

r  * 


'?.**  t  ******  t*t  If  ******t*************1(******  t***tt*****tt*****t***1(  ****** 

c 

SUBROUTINE  CORVEC ( ALPHA, BETA , CORA , CORB) 

INTEGER  I , N 

REAL  ALPHA, BETA, CORA, CORB 

REAL  DERA, DERB, TEMP, RS, WEIBUL, Cl 1 ,C12,C22,D1 ,D2 
REAL  X ( 15  > , Y ( 15  > 

COMMON  N , X  ,  Y 
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10 


Cl  1-0,0 
012  =  0 . 0 
022  =  0 . 0 
D 1  =  0 . 0 
D2  =  0 . 0 


Oil .C12.C22.D1 .  AND  D2  REPRESENT  A  SYMMETRIC 
SYSTEM  OF  2  EQUATIONS  IN  2  UNKNOWNS,.  THE  2 
UNKNOWNS  ARE  CORA  AND  CORD.  HERE  Cll •C12fC22* 
D1 *  AND  D2  ARE  INITIALIZED  TO  0.  IN  THE  DO 
LOOP  THAT  FOLLOWS*  THEIR  VALUES  ARE  COMPUTED, 


DO  10  1  =  1  *N 

RS  =  Y( I >-WEIBULCX< I ) .ALPHA. BET A) 
T  E  M  P  =  X < I >**BETA 
DERA=TEMP*EXP(-ALPHA*TEMP) 

DERB=DERA*ALOG(X( I ) )* ALPHA 

C11  =  C11  +  <DEF'A**2) 
C12=C12+'DERA*DERB) 

022=022+ ( PERB**2 ) 

D1 =D1 + ( DERA*RS ) 

D 2  =  P2+<  BERBERS ) 

CONTINUE 


R  S  =  0  B  S  -  E  X  P 


DERIVATIVE 
TO  ALPHA 
DERIVATIVE 
TO  BETA 
COMPUTE  Oil 
COMPUTE  012 
COMPUTE  022 
COMPUTE  D1 
COMPUTE  D2 


WITH  RESPECT 
WITH  RESPECT 


TEMP* (Cl 1*C22 ) - <C12**2> 


C  0  R  A  =  ( (D1*C22)-(D2*C12> )/TEMP 
CORB=  (  <C11*D2)-(D1*C12)  > /  TEMP 
RETURN 
END 


'  NOW  THE  SYSTEM  IS  SOLVED 
1  USING  CRAMER'S  RULE 
!  TEMP  IS  THE  DETERMINANT 


************##***##*#*##*#*****####*#***###***  ************  i******** 


C  Y 

0  MODLIN  IS  A  SUBROUTINE  THAT  IMPLEMENTS  THE  MODIFICATION  OF  THE  * 

C  LINEARIZATION  METHOD  PROPOSED  BY  H.O.  HARTLEY  IN  HIS  PAPER  ’THE  V 

C  MODIFIED  GAUSS-NEWTON  METHOD  FOR  THE  FITTING  OF  NON-LINEAR  REGRESS  > 
C  ION  FUNCTIONS  BY  LEAST  SQUARES.’  ALL  COMMUNICATION  WITH  THE  PPOC- 
r  DURE  IS  DONE  THROUGH  THE  PARAMETER  LIST?  ALPHA  *  BET A  *  CORA  *  CORE .  + 

C  MODLIN  OPTIMIZES  THE  CORRECTION  VECTORS  COMPUTED  BY  CORVEl  AND  THEN  • 
C  ADDS  THEM  TO  ALPHA  AND  BETA.  WHEN  THE  PROCEDURE  RETURNS*  ALPHA  AND  * 
S  BET A  ARE  THE  NEW  PARAMETER  ESTIMATES.  THE  VALUES  OF  CORA  AND  CORD  * 
C  MAY  HAVE  BEEN  CHANGED  IN  THE  PROCEDURE.  * 

r~  .i 


C  ********  ft  *  *  *  *  *  *  *  *  *  *  *  ********  *  *  *  ***********  *************************** 


SUBROUTINE  MODL IN ( ALPHA .BETA  *  CORA . CORB ) 

REAL  ALPHA. BETA. CORA. CORB 
REAL  QO  *  Q .1.  *  Q2  .  V  *  SSE  *  PENOM 
REAL  TEMP 

!  LET  THETA* ( ALPHA  *  BETA )  BE  THE  CURRENT  PARAMETER  VALUES  AND 
!  DELTA=(CORA*CORB)  BE  THE  CORRECTION  VECTOR.  MODLIN  ESTIMATES 
!  THE  VALUE  OF  V>  =  0  THAT  MINIMIZES  SSE  <  THET A+V*DELT  A  '>  ,  SSE  IS 
i  COMPUTED  AT  THETA  (QO)*  THETA+  .  5*DEL  T  A  (Ql>.  AND  THETA4-  D  E  L  T  A 
’  (Q2).  V  IS  FOUND  SO  THAT  THETA+V* DELTA  IS  THE  VERTEX  OF  THE 

!  PARABOLA  MASSING  THROUGH  QO*  Qt*  AND  Q2. 
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C10  =  SSE<  ALPHA»BETA> 

Gl=SS£(ALPHA+.5*CQRAfBETA+.5*C0RB> 

G2=SSE( ALPHA+CORA.BETA+CORB) 

C 

10  DEN0M=4.0*(Q2+G0-(2.0*01)  ) 

C 

'  IF  DENOM  IS  CLOSE  TO  ZERO  THEN  WE  CAN'T  COMPUTE  V  WITHOUT 
!  PRODUCING  A  DIVIDE-BY-ZERO  OR  FLOATING-POINT-OVERFLOW 
!  ERROR.  IN  THIS  CASE*  ADD  THE  CORRECTION  VECTOR  ASSOCIATED 
1  WITH  THE  MINIMUM  OF  01  AND  02  TO  ALPHA  AND  BETA. 

C 

TF  (ABS'DENOM) .LT. IE-15)  THEN 
IF  (01 .LT ,02)  THEN 
AL.PHA  =  ALPHA+,5*C0RA 
BETA=BETA+.5*C0RB 

ELSE 

ALPHA -ALPHA+ CORA 
B  E  T  A B  F  T  A  F  C  0  R  B 
END  IF 
RETURN 
END  IF 
C 

V .  5  +  <  (00-02) /DENOM) 

TEMP  "SSE(AL.PHA  +  V*CQRA*BETA  +  V*CORB) 

!  IF  V<0  OR  TEMP >00  THEN  THE  COMPUTATION  OF  V  IS  REDONE  WITH 
!  DELTAS. 5*DELTA. 

P 

IF  ( (V.LT.O.O) .OR. (TEMP.GT.QO) )  THEN 
C  0  R  A  = . 5 ♦COR A 
C  0  R  B  =  .  5  *  C  0  R  B 
Q  2-01 

Q1=SSE( ALPHA+.5*C0RA*BETA+.5*C0RB> 

GOTO  10 
END  IF 
r 

ALPHA= ALPHA+ V*CORA 
B  E  T  A  =  B  E  T  A  +  V  *  C  0  R  D 
RETURN 

END 


(7  ******************************************************** ************ 


r  * 

:  GUESS  IS  A  SUBROUTINE  THAT  FINDS  STARTING  VALUES  FOR  ALPHA  AND  ♦ 
r  BETA.  COMMUNICATION  WITH  THE  PROCEDURE  IS  THROUGH  THE  PARAMETER  * 
C  LIST  AND  THE  COMMON  BLOCK,  ALPHA  AND  BETA  ARE  BOTH  OUTPUT  PARA-  * 
C  METERS  f  ACC  IS  AN  INPUT  PARAMETER,  * 
C  * 


c  * it**************************************************. **********  ******* 
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SUBROUTINE  GUESS  <  ALPHA  *  BETA  *  ACC ) 

EXTERNAL  PSSEA  * PSSEB 
INTEGER  I  >  J  >  N  *  E  <  3  )  *  ACC' 

REAL  ALPHA  *  BETA  *  C ( 3  > *D<3) * T 1 , T2 * PSSEA * PSSEB 
REAL  X'l!i'-fY<l?) 

LOGICAL  CONVERGE 
COMMON  N  *  X  *  Y 

1  GIVEN  TWO  DATA  POINTS  IN  THE  EMPIRICAL  DISTRIBUTION*  WE  CAN 
!  SOLVE  FOR  ALPHA  AND  BETA  SO  THAT  THE  WE  I  BULL  MODEL  FIT'.? 

!  THROUGH  THOSE  TWO  POINTS  EXACTLY.  GUESS  CHOOSES  THREE  DA i A 
i  POINTS  AND  FITS  THE  WE  I  BULL  THROUGH  THE  HIRST  TWO  *  THF  LA  S’ 

1  TWO>  AND  THE  FIRST  AND  LAST  POINTS*  THUS  ARRIVING  AT  THREE 
1  DIFFERENT  ESTIMATES  FOR  ALPHA  AND  BETA.  THE  AVERAGE  7  Of  r''l 
'  THREE  ESTIMATES  ARE  USED  AS  STARTING  VALUES  FOR  HPUA  AM" 

!  BETA 

ALPHA =0.0 
BETA=0 . 0 

IF  '.Y(N>  ,  EH,  0.0)  RETURN 

E ( 1 > =2  !  THE  VALUES  OF  E(l>  *  E(2)»  AND  E(3>  DE  TERM T NL 

E  <  2 )  r-8  !  WHICH  THREE  DATA  POINTS  ARE  CHOSEN.  HERE  T  HI 

E  <  3 ) = 1 3  !  2ND  *  8TH  *  AND  13TH  POINTS  ARE  USED, 

DO  20  1=1*3 
D-:i  )=X(E(  I)  > 

IF  ( Y • E< I > > .Ed. 0.0)  THEN 
C<I) =.00001 
ELSE 

C ( I >  =Y(E(  I  >  ) 

END  IF 
CONTINUE 

DO  30  1=3.3 
J  =  MOD  < I  *  3  >  + 1 

T 1 = ALOG  <  ALOG ( 1 . 0~C ( I ) >/ALOG< 1 . 0-C ( J) >  > /ALOG ( IK  I )/IK J  > ) 

DETA  =  BET  A  +  T 1 

ALPHA -"ALPHA* ( -ALOG ( 1 . 0 - C ( T ) > /D ( I >*#T1  ■ 

CONTINUE 

ALPHA= ALPHA/3 . 0 
BETA-BETA/3.0 

1  MOST  OF  THE  TIME  THESE  STARTING  VALUES  WILL  BE  GOOD  ENOUGH  T 
!  BEGIN  THE  NONLINEAR  REGRESSION  PROCEDURE.  SOME  CASES 
'  HOWEVER*  WILL  REQUIRE  EVEN  MORE  ACCURATE  STARTING  VALUES. 

1  THE  FOLLOWING  CODE  OPTIONALLY  IMPROVES  THE  STARTING  VALUES 
'  DEPENDING  ON  THE  VALUE  OF  ACC.  ACC  IS  MERELY  THE 
'  NUMBER  OF  TIMES  THE  LOOP  BELOW  IS  EXECUTED.  THE  LOOP  TP  TEC 
'  TO  OPTIMIZE  ALPHA  FOR  A  FIXED  BETA*  AND  THEN  OPTIMIZES  BETA 
!  FOR  A  FIXED  ALF'HA  . 
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40 

so 


r 


IF  'ACC.LE.O)  RETURN 

DO  40  I"1 , ACC 
T 1  =ALPHA 
T  2  B  E  T  A 

CALL  SECANT ( ALPHA , BETA , ALPHA ,PSSE A » CON VERGE) 

IF  (CONVERGE)  CALL  SECANT ( ALPHA , BET A , BETA , PSSEB .CONVERGE ) 
IF  ( .NOT. (CONVERGE) )  GOTO  50 
CONTINUE 
RETURN 
ALF'HA  =  T  1 
B  E  T  A  =  T  2 
RETURN 
END 


C'¥V¥¥Y¥¥'¥'+#*>f*  +  *)¥)¥*:¥*)¥)¥*)¥)¥>¥)¥)¥)¥**#>¥*i¥¥*>¥>¥>¥*****)¥**1¥  ************  *Y'¥  *#¥■**'¥ 

* 

r  FAIL  IS  AN  OUTPUT  ROUTINE  CALLED  WHEN  A  DISTRIBUTION  FAILS  TO  CON-  * 

C  VERGE  AFTER  50  ITERATIONS.  THE  VALUES  OF  SEVERAL  VARIABLES  ARE  * 

C  WRITTEN  TO  THE  OUTPUT  FILE.  COMMUNICATION  WITH  THE  PROCEDURE  IE  * 

0  DONE  THROUGH  THE  PARAMETER  LIST  AND  THE  COMMON  BLOCK .  ALL  PARA-  ¥ 

C  METERS  ARE  INPUT  PARAMETERS.  ¥ 

C  1 

C  ¥  '¥  ¥  *  ¥  ¥  ¥  :¥**¥¥'¥¥:¥  ¥  ♦  *  *  *  *  ¥  *  t  *  ¥  % %%%%%%%%%%%*■*}%%%  *•*,%%**  % :*  #  **#**)¥*  *  ¥  ¥  *  ¥  ¥  ¥  ¥  ¥  ¥ 


SUBROUTINE  FAIL ( SID . MONTH. HOUR . STARTA. STARTS . ALPHA . BETA . NRSS > 
INTEGER  SIDf MONTH fHOURfNf I 

REAL  START A ,  ST ARID , ALPHA . BET A , NRSS  , TEMP .X<15>.YU5>. SEE 
COMMON  N.X.Y 


TEMP "SSE'ETAPTA- STARTS) 

W  RITE  >’  6  ,  100  '• 

WRITE < 5 ■ 200  >  SID . MONTH  » HOUR 

Wr  1  T  E '  L  .  300  > 

WRITE ( A  . 400 )  STARTAf STARTS . TEMP 
WRITE (6.500)  ALPHA. BET A . NRSS 
WRITE (6. 600) 

WRITE (6 .700) 

DO  10  1=1. N 

WRITE' A . S  0  0  >  X(I)  ?Y(I> 

10  CONTINUE1 

RETURN 


1.00  FORMAT  <  /  •'/  .  T32 .25  <  '  ¥  '  )  ./// ) 

:”>0  FORMAT  (  'rlX*  'ATTENTION?  STATION  '.12.'.  MONTH  -  .12.  '  f  HOUR  '• 
‘T  11.'  FAILED  TO' ./.  IX.  'CONVERGE  AFTER  50  ITERATIONS.') 

300  FORMAT'/. IX. 'A  VARIABLE  DUMP  FOLLOWS:'./) 

1 00  FORMAT' IX F ' ST  ART A*  '  >015.7, '  STARTS* ' .G15.7- 
t  '  SSE(STARTA, STARTS)*' ,G15. 7) 

COO  FORMAT( IX, ' ALPHA= ' ,G15. 7, '  BET A* ' , G 1 5 . 7 , 

*  '  SSE(ALPHA,BETA)=' ,G15.7) 

fOO  FORMAT'/, 5X,  '  ENDPTS  '  ,14X,  '  OBCUMFF: '  ) 

’00  FORMAT  (  '  +  '  ,  4X  ,  ' _ '  ,  14X,  ' _  >  /> 

SOO  FORMAT' 1X.G15. 7.5X.G15.7) 

F  N  0 
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r*:tV******r*  ******♦***#*******♦************#**+***#****  *******  +  ****♦** 


'.:  SUCCESS  IS  AN  OUTPUT  ROUTINE  CALLED  WHEN  A  DISTRIBUTION  CONVERGE  •  t 
C  SUCCESSFULLY  WITHIN  50  ITERATIONS.  SUMMARY  STATISTICS  OF  T  HF.  K 
C  REGRESSION  ARE  WRITTEN  TO  THE  OUTPUT  FILE.  COMMUNICATION  WITH  i 
C  THE  PROCEDURE  IS  THROUGH  THE  PARAMETER  LIST  AND  COMMON  BLOCK  * 
C  ALL  PARAMETERS  ARE  INPUT  PARAMETERS.  * 


c  ¥ 

r#*  ******* ***** *  * **** ****** **** ******* ***** ************ *#* ****** *.**  * 

SUBROUTINE  SUCCESS ( S I D > MONTH > HOUR » ALPHA > BET A . NRSS > COUNT > 

INTEGER  SI  Dr  MONTH > HOUR > I >N .COUNT 

REAL  AL.  PHA  «  BETA  >  NRSS  >  RMS  >  X  (  1 5  >  >  Y(  15)  >  EX  >  RS  »  WE  I  BUI. 

COMMON  NfXfY 

RMS=SGRT ( NRSS/FLOAT ( N ) ) 

WRITE',  A?  100) 

WRITE <6 >200> 

WRITE <  6  >  300  > 

WRITE ( 6  >  A 0 0 )  SID >  MONTH  > HOUR  >  ALPHA  >  BETA  > RMS >  COUNT 
WRITE (5 >500' 

WRITE! 6>600> 

DO  10  1  =  1  >M 

EX=WEIBUL (X  ( I ) > ALPHA > BETA) 

RS-Y( I ; -EX 

WR I TE ( 6  v 700 )  X ( I ) >  Y •  I > >  EX  > RS 
CONTINUE 
RETURN 

FORMAT < /// »  T31 »  35  < ’*  ')>///) 

FORMAT ( 1  OX  *  ' STATION  ID" >  5  X  > 'MONTH' >5X>  'HOUR ' • 1  OX r ' AL PH A  '  ' 1 5 » 

*  'BETA' >17X> 'RMS' >11X> '*  OF  ITERATIONS') 

FORMAT  (  '  +  '  >  9  X  >  ' _ __'>5X>' _ '  »5X>  ' _ '  >  10X>  ' .  '  >15X» 

t  ' _ '  >17X>  ' _ '  >  11X>  '_  __  '  >//> 

FORMAT  <  12X  >  I  &  >  2  ( 8X  >  12  )  >  6X  >  G14 . 7  >  4X  >  G  3  4  .  7  >  4X  >  G 1,  4 . 7  >  i  IX  >  1 3  >  V  '• 

FORMAT <T38> 'ENDPTS ' >5X>  ' OBCUMFR ' > 10X » 'EXCUMFP ' >  12X >  'RESIDUAL  '  > 

FORMAT  <  '4  '  >T38>  ' _ '  »5X>  ' _ '>10Y,' _ '»12X-  .  . 

*,  /  /  > 

FORMAT ( T37  * F7 . 4  >  AX  >  F5 . 3  >  7X  >  G15 . 7 >  5X  »  G15 . 7 > 

LND 

t-*********  **!»;*  *******  ******************************************  ******** 


C  Y 
C  SECANT  IS  A  SUBROUTINE  THAT  USES  THE  SECANT  METHOD  OF  ROOT  SOLVING  ♦ 
C  TO  FIND  THE  ROOT  OF  PSSEA  HOLDING  BETA  CONSTANT.  OP  TO  FIND  THE.  t 
C  C,OOT  OF  PSSEB  HOLDING  ALPHA  CONSTANT.  (THIS  MEANS  IT  WILL  FIND  THE  * 
0  BEST  ALPHA  FOR  A  GIVEN  BETA.  OR  THE  BEST  BETA  FOR  A  GIVEN  ALPHA.)  ■> 
r  COMMUNICATION  WITH  THE  PROCEDURE  OCCURS  THROUGH  THE  PARAMETER  LIST  * 
C  AND  THE  COMMON  BLOCK.  IN  THE  PARAMETER  LIST.  ALPHA  AND  BETA  ARE  * 
"  THE  CURRENT  VALUES  OF  THE  MODEL  PARAMETERS.  FARM  IS  THE  VARIABLE  > 

.  fr:Nr.  OPTIMIZED  (EITHER  ALPHA  OR  BETA )  >  AND  PDER  IS  THE  PARTIAL  * 
C.  DERIVATIVE  FUNCTION  (EITHER  PSSEA  OR  PSSEB).  CONVERGE  IS  A  LOGICAL  * 
C  VARIABLE  INDICATING  WHETHER  THE  SECANT  METHOD  CONVERGED  ON  A  ROOT.  * 
0  *  V  SOLVE  -OR  THF.  ROOT  OF  PSSEA  >  SET  PARM-ALPHA  AND  PDER=PSSEA,  THE  % 


r  OPTIMIZED  VALUE  OF  ALPHA  WILL  BE  RETURNED.  TO  SOLVE  FOR  THF  ROOT  OF* 
f  PSSEB >  SET  p A R M  ~  B  E  T  A  AND  PPER=PSSEB.  THE  OPTIMIZED  VALUE  OF  BETA  * 
r  WILL  RF  RETURNED.  » 

* 

*******  ******************************************************* 


1  o 
r 

100 

200 

300 

ROO 

500 

r, 
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SUBROUTINE  SEC  ANT  (  ALPHA p BETA » FARM » PI'ER  p  CONVERGE  ) 

REAL  ALPHA p BETA p PARM p PDER p X < 15) pY( 15) pTI pT2pT3p DELTA 
INTEGER  Nr  I 
LOGICAL  CONVERGE 
COMMON  N  p  X  p  Y 

r 

T2=PDER< ALPHA p BETA) 

DELTA® • 00 1 
P ARM® PARM+ DELTA 
CONVERGE® . TRUE . 

r 

DO  10  1=1 p 15 

Tl® PDER (ALPHA. BETA) 

T  3  ~  T 1  -  T  2 

IF  C  A  B  S  <  T  3 ) >LE. IE-15)  GOTO  15 
DELTA- ( - T 1 ♦DEL TA ) / ( T3 ) 

FARM®  F'APM+DELTA 

IF  (ABS( DELTA) .LE.1E-7)  GOTO  20 
T2=T1 

i.0  CONTINUE 

15  CONVERGE® . FALSE  . 

20  RETURN 

END 
r 

C  f  ♦  ♦  ♦♦♦♦♦♦♦♦  *'*  ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦•*♦♦♦♦♦♦♦♦♦♦♦♦4  ♦♦*» 

I"  * 

C  PSSEA  IS  A  REAL  FUNCTION  THAT  COMPUTES  THE  PARTIAL  DERIVATIVE  OF  SSL* 
•:  HUH  RESPECT  TO  ALPHA.  COMMUNICATION  WITH  THE  PROCEDURE  OCCURS  # 

THROUGH  THE  FUNCTION  NAME  p  COMMON  BLOCKp  AND  PARAMETER  LIST,  ALPHA  * 
C  AND  BETA  ARE  INPUT  PARAMETERS?  THE  VALUE  OF  THE  DERIVATIVE  IS  "f 

■S  RETURNED  THROUGH  THE  FUNCTION  NAME,  * 

A 

r*******  ♦  ♦♦♦♦♦  ♦  .?***.#  ♦♦*♦♦♦♦♦*♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦!((♦  ♦♦♦♦♦♦♦♦♦♦♦♦»«*%**¥*$ 

REAL  FUNCTION  PSSEA ( ALPHA p BET A ) 

REAL  ALPHA  r  BETA  p X  < 1 5 ) . Y  < 1 5  )  p  T 1 p T? 

INTEGER  Nr  I 
COMMON  NrXrY 

r 

PSSEA® 0,0 
DO  10  1=1. 1’ 

T 1 = X (  I  ) *tPET A 
IS " EXP ( -ALPHA+T 1 ) 

PSSEA=PSSEA+  <  Y  I)-l  ,  0  +  T2 )  ♦  ( -T 1  *T2  > 

10  continue 

r 

RETURN 
E  M  D 
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C  *  *  ?  *  *  *  *■  *  *  *  *  *  »  ¥  *  *  *  *  *  * *  *  ¥  *  *  *  *  *  *  *  *  *  *  *  *  Dr  *  *  Hr  *  If  %  *  *  *  *  *  t  *  *  *  *  t  *  t  ¥  ¥  *  *  *  *  <f  a  *  *  *  M  ♦  * 


c  > 
C  PSSEB  IS  A  REAL  FUNCTION  THAT  COMPUTES  THE  PARTIAL  PEP T VAT ! UP  OF  * 
0  SSE  WITH  RESPECT  TO  BETA.  COMMUNICATION  WITH  THL  PROCEDURE  OCCURS  * 
r  THROUGH  THE  FUNCTION  NAME  t  COMMON  BLOCK*  AND  PARAMETER  LIST  ALPHA  * 
C  AND  FETA  ARE  INPUT  PARAMETERS?  THE  VALUE  OF  THE  DERIVATIVE  I '  V 
C  RETURNED  THROUGH  THE  FUNCTION  NAME.  1 


C*********#*****#**#***#******************#***?**?*  *♦♦**»♦*♦♦**♦  *'**¥»** 

c 

REAL  FUNCTION  PSSEB < ALPHA . BETA > 

REAL  ALPHA  »  BETA »  X ( 15  > *Y(15)*TlrT2*T3fT4 
INTEGER  N » I 
COMMON  N»XfY 

r* 

PSSEB-0.0 
DO  10  I-lrN 
TL-X( I > ¥ ¥  BETA 
T2-EXP i - ALPHAKT 1 ) 
v3- -ALPHA* AL.OG  (  X  ( I )  >*T1*T2 
T4=Y( I >  - ! .0+12 
F'SSEB-PSSEB  +  T3*T4 

:i  0  CONTINUE 

r 

RETURN 

END 
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NONLINEAR  REGRESSION  OF  THE  UEIBULL  MODEL  ON  VISIBILITY  DATA 
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4.  Some  Remarks  About  the  Program 


The  program  given  in  this  report  was  written  to  fit  the  Weibull 
distribution  to  visibility  data.  The  visibility  data  was  that 
contained  in  the  "Revised  Uniform  Summary  of  Weather  Observations" 
(RUSSWO's)  prepared  by  the  Data  Processing  Division  of  the  Air  Weather 
Service.  Some  changes  must  be  made  to  fit  another  distribution  to 
another  variable.  The  program  is  made  up  of  a  series  of  subroutines 
and  functions  so  that  these  may  be  altered  to  fit  the  users  need  without 
changing  the  flow  of  the  program. 

The  function  WEIBUL  must  be  replaced  by  the  desired  cumulative 
distribution  function  (CDF).  Also,  the  name  of  the  function,  WEIBUL, 
must  be  changed  throughout  the  program  if  the  name  of  the  function  is 
changed.  The  subroutine  GUESS  which  gives  the  initial  values  of  the 
parameters  must  be  changed  to  correspond  to  the  distribution  being 
used.  The  same  basic  idea  can  be  used.  However,  it  may  take  more 
programming  for  other  distributions  particularly  if  the  distribution 
is  not  in  closed  form.  The  subroutine  CORVEC  must  be  changed  where 
the  partial  derivatives  DERA  and  DERB  occur.  Another  change  is  required 
in  the  functions  PSSEA  and  PSSEB  which  are  functions  which  take  partials 
with  respect  to  each  of  the  parameters  of  the  sum  of  squared  errors. 

The  appropriate  partials  must  replace  those  of  the  Weibull  distribution 
in  each  of  these  functions. 

The  program  was  designed  to  run  in  a  batch  environment,  and  the  rules 
of  standard  FORTRAN  -  IV  were  adhered  to  as  closely  as  possible.  The  only- 
departures  from  FORTRAN  -  IV  are  the  use  of  the  IF-THEN-ELSE-ENDIF 
structure  commonly  found  in  FORTRAN  77,  and  the  use  of  the  exclamation  mark 
to  permit  comments  on  the  same  line  as  code.  Both  deviations  were  made 
merely  for  clarity's  sake  in  the  program. 


5.  Sample  Output 

A  sample  output  is  shown  in  page  14. 
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